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1.  Introduction 

Hyperbolic  partial  differential  equations  typically  contain  regions  of  greater  activity  such  as 
shocks  or  wave  fronts,  and  these  features  are  dynamic.  This  behavior,  combined  with  the  difficulty 
of  providing  a  global  uniform  grid  fine  enough  to  resolve  the  smallest  important  length  scale,  has 
stimulated  much  work  on  adaptive  grid  techniques  for  hyperbolic  problems.  There  are  two  major 
approaches.  The  first  is  to  use  an  irregular  grid,  where  the  choice  of  grid  points  depends  on  some 
function  of  the  solution.  The  grid  may  be  continuous  in  time.  The  second  is  to  use  locally  uniform 
grids,  and  to  replace  these  grids  frequently.  In  this  case  the  grids  are  not  continuous  in  time. 

This  second  method  is  called  Local  Uniform  Mesh  Refinement  (LUMR)  and  is  a  powerful 
technique  for  the  solution  of  partial  differential  equations.  It  consists  of  adaptively  refining  a 
uniform  coarse  grid  by  overlaying  it  with  uniformly  refined  subgrids.  The  uniformity  is  important 
in  insuring  that  the  overhead  of  managing  the  refinement  is  kept  small,  and  that  the  refinement 
does  not  defeat  any  vectorizability  which  may  be  present  in  the  discretization.  This  method  has 
been  shown  to  be  efficient  for  a  variety  of  one  and  two  dimensional  time  dependent  problems  [1, 
4,  7].  We  discuss  how  to  combine  LUMR  and  moving  grids  in  a  way  which  substantially  improves 
the  performance  of  each  of  these  algorithms  at  very  little  additional  cost. 

In  Section  2,  we  motivate  this  extension.  In  Section  3,  we  give  a  detailed  description  of  the 
algorithm.  In  Section  4,  we  show  computational  experiments  of  this  method,  and  compare  it  with 
uniform  grid  and  regular  LUMR  calculations. 

2.  Motivation 

In  a  numerical  approximation  to  a  hyperbolic  PDE,  the  truncation  error  typically  has  leading 

form 

r=C7(t)*(A'  +  /fc')  (2.1) 

where  C(t)  is  independent  of  A  and  k.  The  usual  approach  in  reducing  the  local  error  is  to  reduce 
A  and  k.  However,  since  p  is  usually  small  ((8]  argues  that  p  =*  4  is  optimal  for  u »  =  «,),  this  is 
only  moderately  efficient.  In  particular,  if  the  problem  has  n  space  dimensions,  the  work  goes  as 
(assuming  k/h  *  A  is  constant)  A“("+‘l.  Thus,  to  reduce  the  global  error  by  a  factor  of  2  by  using 
mesh  refinement  will  require  a  decrease  in  h  and  k  by  2l/p  and  an  increase  in  the  work  of  2^n+U/r. 
Thus,  the  small  p  common  in  the  solution  of  hyperbolic  problems  makes  mesh  refinement  relatively 
inefficient  at  reducing  the  error  in  the  computation. 

What  alternatives  are  there?  It  is  also  possible  to  reduce  the  local  error  by  reducing  C(t)  in 
(2.1).  We  may  be  able  to  do  this  by  a  change  of  variables  to  a  moving  coordinate  system,  as  C(t) 


depends  on  the  problem  being  solved.  For  example,  consider  the  problem 

+  ati,  =  0,  (2.2) 

approximated  with  the  leap-frog  scheme.  The  principle  truncation  error  is 

T  m  |(**u,«  +  aA8u,„)  (2.3) 

We  can  use  (2.2)  to  eliminate  u«{  from  (2.3)  to  get 

T  -  gOU,„(-a8*s  +  As)  (2.4) 

so  that  C(t)  — ►  0  as  a  — *  0.  It  is  easy  to  pick  a  coordinate  transformation  (x,  t)  — »  (x,  t)  which 
takes  (2.2)  to 

u,-  +  0  •  us  =  0. 

Moreover,  leap-frog  is  exact  for  this  problem. 

The  important  thing  to  note  in  (2.4)  is  that  reducing  a  by  2  reduces  the  truncation  error  by 
(roughly)  2.  Thus,  a  significant  savings  in  the  effort  required  to  reach  a  given  error  tolerance  can 
be  realized  by  finding  an  approximate  transformation. 

This  is  not  a  new  approach;  part  of  the  success  of  the  moving  finite  element  method  [9,  10] 
can  be  traced  to  this,  as  can  the  method  of  Davis  and  Flaherty  [5].  What  is  new  here  is  a  way 
to  reduce  C(t)  in  some  cases  at  little  additional  cost  in  a  mesh  refinement  code.  We  will  not  try 
to  drive  C{t)  to  zero  (though  that  is  possible  in  exceptional  cases),  rather,  we  will  apply  uniform 
velocity  transformations  to  each  refined  grid  in  an  attempt  to  reduce  C(t)  at  low  cost. 

We  also  argue  below  that  this  is  an  effective  way  of  reducing  numerical  dispersion  in  discon¬ 
tinuous  or  nearly  discontinuous  solutions.  The  best  way  to  see  this  is  to  look  at  graphs  of  the 
numerical  dispersion  of  a  sample  scheme  on  a  model  problem.  We  consider  the  leap-frog  method 
applied  to  u<  +  au,  =  0  for  0  <  a  <  1.  The  model  is  U{  +  u,  =  0  where  a  here  is  the  result  of 
shifting  to  a  moving  coordinate  system.  We  let  k/h  =  |.  This  gives  the  dispersion  relation 

wA  =  2sin~l  ^asin(A^) 

where  the  basic  mode  is  exp(i(wt  -  £z)).  However,  this  is  not  quite  what  we  want  to  compare. 
Because  we  are  using  a  moving  coordinate  system,  we  actually  have  (in  the  original,  stationary 
frame) 

wA  m  2sin_I  f  iasinfA^A  +  (1  -  a)A{ 


Figure  1:  Dispersion  relation  for  model  problem  for  various 
values  of  a.  Parasitic  waves  have  been  omitted. 

Graphs  of  this  for  various  values  of  a  are  shown  in  Figure  I. 

The  importance  of  this  is  that  even  for  a  —  the  numerical  dispersion  is  significantly  re* 
duced.  Thus,  even  an  approximate  moving  transformation  can  be  very  beneficial.  Also  note  that 
this  reduction  in  numerical  dispersion  extends  right  to  the  highest  frequencies,  so  this  applies  to 
discontinuous  solutions  as  well. 

3.  Algorithm 

In  this  section  we  discuss  the  algorithm  for  moving  LUMR  (MLUMR)  in  two  space  dimensions. 
We  first  discuss  the  grid  structure,  then  the  regridding  strategy,  and  finally  the  algorithm  for 
advancing  the  solution  one  time  step.  Our  algorithm  is  more  in  the  flavor  of  [7]  than  [1]  in  that 
we  don’t  allow  the  grids  to  be  arbitrarily  oriented. 

The  algorithm  is  fully  recursive,  so  multiple  levels  of  refinement  are  easily  accommodated. 

3.1.  Grid  Structure 

The  mesh  refinement  algorithm  produces  a  series  of  level »  of  grids,  G(.  The  initial,  user 
specified,  level  (corresponding  to  the  spatial  domain  of  the  problem)  is  denoted  Go-  All  levels  are 
made  up  of  rectangular  pieces,  called  grids,  denoted  Gu  for  the  Ith  grid  in  level  l.  Further,  each  of 
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Figure  2:  Sample  grids.  One  level  of  refinement  is  present. 

these  pieces  is  lined  up  with  the  usual  cartesian  directions.  Finally,  grids  must  be  strictly  nested: 
if  a  grid  Gu  lies  in  Gt~ij,  then  all  of  Gu  lies  in  Cj_ ij  and  none  of  Gu  lies  in  any  other  grid  at 
level  l  -  1.  The  grid  Gj_tj  is  called  the  parent  of  Gu,  and  Gu  a  child  of  Gi-ij.  A  sample  grid 
structure  is  shown  in  Figure  11. 

Note  that  the  grids  Gu  and  GI2  are  two  separate  grids  rather  than  a  single  grid  lying  partly 
in  Got  and  partly  in  G<».  This  restriction  significantly  simplifies  the  data  structures  involved  in 
managing  the  grids  without  imposing  a  restriction  on  the  kinds  of  refinement. 

Each  grid  G<r  contains  a  uniform  mesh  of  points  (2,-,yy),  i  —  l,...,n,  j  =»  l,...,m,  and 
*i  =  *o  +  (<  -  l)h|,  yy  *  y0  +  (j  -  1)A/,  where  (x0,  y0)  give  the  coordinates  of  the  lower  left  corner 
of  G(r.  The  step  sizes  in  space  are  Aj  and  the  step  sizes  in  time  are  k(.  If  the  boundaries  of 
two  grids  abut,  each  contains  the  mesh  points  along  the  common  boundary.  Thus,  points  along  a 
common  boundary  occur  in  two  different  places  in  memory;  the  algorithm  will  show  how  they  are 
maintained  with  a  single  value.  To  simplify  the  programming,  grids  may  only  abut  along  rows. 

It  is  important  to  note  that  the  grids  of  one  level  overlay  the  grids  of  the  coarser  levels,  rather 
than  being  patched  into  the  coarser  grids.  This  maintains  the  uniform  nature  of  the  grids,  at  a 
small  cost  in  additional  grid  points.  As  the  same  point  in  the  domain  may  lie  in  several  grids, 
the  algorithm  must  take  steps  to  insure  a  single  valued  solution.  Because  the  grids  are  regularly 
oriented,  this  can  be  done  in  an  efficient,  vectorizable  fashion. 

Each  grid  is  also  assigned  a  velocity  ff.  The  entire  grid  will  be  moved  uniformly  at  this  velocity, 
with  v  changing  as  often  as  every  time  step.  Because  the  grids  move,  they  can  run  into  each  other 
and  escape  from  their  parent. 

In  the  first  case,  there  are  two  possibilities:  one  grid  is  overtaking  the  other,  or  two  grids 
run  into  each  other.  In  the  case  of  one  grid  overtaking  another,  we  allow  the  grids  to  overlap. 
Nothing  else  need  be  done,  as  this  will  still  provide  a  good  local  match  to  the  solution  velocity. 
Alternately,  we  could  force  the  grids  to  abut,  and  reset  their  velocities  so  the  area-average  of  their 
respective  velocities.  In  the  ease  where  two  grids  collide,  this  means  that  there  is  no  local  velocity 


4 


(a)  (b) 

Figure  S:  Colliding  grids  as  in  (a)  are  handled  by  changing 
to  a  stationary  grid,  with  additional  levels  of  refinement  as 
necessary,  as  in  (b). 

transformation,  and  the  error  must  be  reduced  by  reducing  the  step  sites  h  and  k.  Thu  can  be 
handled  naturally  by  MLUMR  by  regridding  and  generating  an  additional  level  of  refinement  as 
necessary  (cf.  Figure  3). 

In  the  second  case  of  grids  escaping  from  their  parent  into  another  grid  as  in  Figure  4,  we 
force  the  grid  to  stop  at  the  parent  boundary,  and  regrid  to  create  a  grid  in  the  adjacent  parent. 
Because  the  grid  can’t  move  (relative  to  its  parent),  it  may  be  necessary  to  refine  it  as  well;  again, 
this  is  handled  in  a  natural  fashion. 

Another  approach  to  this  problem  is  to  sacrifice  strict  nesting  by  parent  and  go  to  the  looser 
level  nesting  discussed  in  [2].  This  complicates  the  data  structures.  Because  of  the  grid  organi¬ 
zation,  only  grids  at  the  first  level  of  refinement  are  likely  to  escape  from  their  parent  (a  grid  in 
the  coarsest  grid).  The  coarsest  level  will  normally  have  very  few  components — only  one  in  many 
cases.  Grids  at  higher  levels  of  refinement  are  likely  to  move  with  the  same  speed  as  their  parent. 
Thus  grids  will  rarely  escape  from  their  parents,  and  the  additional  cost  of  regridding  when  that 
happens  is  compensated  for  by  the  simplicity  of  handling  the  other  grid-parent  interactions. 

3.2.  Regridding 

After  some  number  of  time  steps  on  a  level  l,  it  is  necessary  to  check  to  see  if  the  refined 
grids  at  the  finer  levels  need  to  be  replaced  with  new  refined  grids.  This  is  done  in  three  steps. 
First,  points  in  the  grids  at  level  /  are  flagged  if  it  is  necessary  to  refine  about  that  point.  This 
may  be  done  by  estimating  the  truncation  error  at  that  point  or  by  some  ad-hoe  method  based 
on  user  knowledge  of  the  behavior  of  the  solution.  Second,  these  flagged  points  are  surrounded 
by  a  buffer  of  marked  points.  This  buffer  is  needed  to  isolate  the  interior  of  the  refined  grid  from 
the  coarse  grid,  which  is  necessary  given  our  strategy  of  overlaying  the  refined  grids.  Further,  this 
buffer  allows  us  to  take  more  steps  between  regriddings  by  giving  a  margin  for  the  region  of  higher 
error  to  move  around  in.  Finally,  the  marked  and  Sagged  points  are  scanned  to  determine  the  new 
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(a)  (b) 


Figora  4:  Grids  leaving  their  parent  as  in  (a)  are  handled 
by  using  stationary  grids  in  the  original  parent  and  the  new 
parent.  An  additional  level  of  refinement,  as  in  (b),  may  be 
necessary. 


Figure  5:  Regridding  process.  •  denote  mesh  points,  •  denote 
flagged  points,  x  denote  buffer  points.  The  resulting  refined 
grids,  G(i  through  Gu  are  outlined. 

refined  grids,  and  those  grids  are  initialized  with  values  from  the  previous  refined  grids  (where  they 
overlapped)  and  the  parent  grid. 

The  process  of  adding  the  buffers  can  be  done  by  bit-wise  logical  operations  of  left  and  right 
shifts  and  or’s.  This  helps  keep  the  cost  of  the  regridding  down.  The  regridding  process  in  shown 
in  Figure  5. 

13.  rinding  grid  velocities 

This  is  either  the  hardest  or  the  easiest  part,  depending  on  what  you  plan  to  do.  In  the  moving 
finite  element  (MFE)  approach,  the  nodes  move  at  a  velocity  which  is  automatically  calculated 
to  minimize  a  penalty  function  made  up  of  the  residual  and  some  terms  to  keep  the  mesh  well 
behaved.  In  our  approach,  we  ask  the  user  to  specify  the  velocity  of  a  grid,  given  the  values  on  the 
grid  and  its  location.  For  many  problems  this  is  quite  reasonable;  detecting  near-discontinuities 
and  computing  an  appropriate  approximate  velocity  is  often  quite  simple.  In  fluids  problems,  the 
appropriate  velocity  is  often  the  fluid  speed,  and  is  immediately  available.  It  should  be  emphasized 
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(*)  (b) 

Figure  0:  Grids  over  which  the  velocity  varies  greatly  as  in 
(a)  are  broken  into  several  grids.  In  (b),  the  grid  has  been 
broken  into  three  overlapping  parts,  with  the  center  part  (in 
thicker  lines)  remaining  stationary. 

that  because  we  are  not  trying  to  get  an  exact  value  for  the  velocity,  we  can  be  much  more  cavalier 
about  choosing  the  velocity.  For  a  completely  automatic  method,  one  similar  to  that  in  [9]  could 
be  used. 

S.4.  Integrating  a  level 

Each  level  consists  of  a  number  of  possibly  intersecting  grids.  Here  we  describe  how  a  single 
level  is  integrated;  in  the  next  section,  we  put  everything  together. 

The  boundaries  between  the  a  refined  grid  and  its  parent  are  determined  by  using  linear 
interpolation  from  the  parent  grid  values.  This  has  been  shown  to  be  stable  [3]  for  stationary 
grids. 

If  the  grids  abut,  then  the  boundary  between  them  is  handled  in  a  special  way.  The  points 
along  the  common  boundary  are  duplicated  in  each  grid.  This  is  done  to  maintain  the  uniform 
structure  of  the  grids.  When  integrating  the  top  grid,  the  common  boundary  is  also  integrated, 
using  data  from  both  grids.  Then  the  values  computed  on  the  common  boundary  are  copied  from 
the  upper  grid  into  the  lower  grid. 

Overlapping  grids  require  special  care.  In  the  case  where  the  refinement  is  1  (no  refinement), 
overlapping  grids  can  be  handled  efficiently  by  just  injecting  the  values  on  the  ‘refined*  grids  onto 
the  parent  grid,  and  then  using  the  usual  boundary  procedure  for  fine/coarse  grid  interfaces  to 
communicate  the  values  back  to  the  refined  grids. 

If  the  refinement  is  greater  than  1,  then  the  boundary  values  must  be  computed  by  interpo¬ 
lating  from  the  overlapping  “partner*. 

As  an  example,  we  will  describe  how  one  time  step  is  taken  on  a  coarse  grid,  where  there  is 
only  one  level  of  refinement  and  no  regridding  takes  place  during  this  time  step.  First,  the  coarse 
grid  is  integrated  from  time  t  to  t  +  Jfeo-  Next,  the  refined  grids  are  integrated  from  time  t  to  t  +  kt, 
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procedure  step(  lev  ) 

1.  Integrate  level  lev 

2.  Check  for  time  to  regrid  the  next  level.  If  so  regrid  the  next  level 

3.  Update  grid  velocities  if  enough  steps  have  passed.  Must  check  this  after  the  regridding  check, 
since  we  don’t  want  to  reset  velocities  and  immediately  regrid. 

4.  Integrate  the  finer  levels,  if  any 
do  i*l.  ref inement (lev) 

call  step(  lev  ♦  1  ) 

5.  Having  integrated  the  finer  level,  we  can  update  this  level  from  the  next  finer  level 

Algorithm  1:  Moving  Mesh  Refinement  algorithm  for  ad¬ 
vancing  the  solution  one  time  step. 

then  from  t+k\  to  t+2ki,  and  so  on  until  they  reach  t+fa.  The  boundary  values  at  the  fine/coarse 
grid  interfaces  are  taken  from  interpolation  of  the  coarse  grid  values  at  times  t  and  t  +  ko-  Finally, 
the  values  on  the  fine  grid  at  time  t  +  ko  are  injected  into  the  coarse  grid  to  update  the  solution 
there. 

3.S.  Integration 

Now  that  we  have  the  structure  of  the  grids  and  the  regridding  algorithm,  we  can  describe 
the  algorithm  for  integrating  a  level  one  step  in  time.  This  algorithm  is  recursive. 

The  algorithm  described  in  Algorithm  1  is  fairly  straightforward.  Step  5  is  the  only  one  we 
have  not  talked  about.  This  step  insures  that  the  solution  on  all  the  grids  is  single-valued,  and, 
more  importantly,  that  those  areas  of  the  coarse  grid  which  are  refined  take  their  solution  from 
the  refined  grids.  If  this  were  not  done  at  every  step,  there  would  be  a  danger  that  the  inaccurate 
solution  on  the  coarse  grid  in  this  area  could  escape  and  contaminate  the  solution  everywhere. 
This  contamination  is  prevented  by  a  combination  of  finite  domain  of  dependence  of  the  differential 
equation  and  the  difference  approximation,  the  buffers  surrounding  the  areas  of  large  error,  and 
the  injection  of  the  better  solution  computed  on  the  fine  grid  into  the  coarser  grids. 

A  more  subtle  point  is  the  order  of  the  steps  in  the  algorithm.  When  the  routine  to  regrid 
is  called,  the  parent  level  has  been  integrated  one  step  beyond  the  step  where  the  regridding  is  to 
take  place.  This  lets  us  place  grids  where  they  will  be  needed,  rather  than  where  they  had  been 
needed.  Also  note  that  when  the  regridding  routine  is  called,  all  grids  at  finer  levels  no  longer  have 
a  legitimate  parent  path  (since  they  now  point  to  non-existent  parents).  This  is  acceptable,  since 


those  grids  will  also  be  regrided  before  they  will  be  used  to  update  any  parent  grid,  and  they  will 
not  be  integrated  (their  only  use  is  to  initialize  the  new  grids) 

Because  handling  overlapping  grids  is  complicated  and  inefficient,  particularly  if  every  grid 
had  a  different  velocity,  we  may  want  to  force  adjacent  grids  which  have  nearly  the  same  velocities 
to  have  the  same  (area  averaged)  velocity.  This  can  significantly  improve  the  performance  of  the 
algorithm.  We  used  this  approach  in  our  program. 

4.  Experiments 

These  experiments  are  all  in  two  space  dimensions  and  for  scalar  problems.  Each  contain 
discontinuities  in  either  the  solution  or  its  derivatives.  Each  was  run  on  a  FPS-164  running  SJE 
and  using  the  E  Release  of  the  FORTRAN  compiler  at  optimization  level  3.  Each  problem 
computed  using  a  uniform  coarse  and  a  uniform  fine  grid,  and  both  LUMR  and  MLUMR.  In  eac 
of  the  mesh  refinement  calculations,  the  level  zero  grid  (the  coarsest  grid)  was  the  same  as  tl 
coarse-only  calculation.  The  width  of  the  finest  mesh  in  each  calculation  is  denoted  by  hj.  N«. 
more  than  1  level  of  refinement  was  used.  Errors  in  the  solutions  were  measured  on  the  coarse 
grid.  Three  different  measures  were  used.  The  ^ -error  is  defined  as  (n-1  e*)1^2.  the  fperror 

is  defined  as  n~1  je;|,  and  the  4»-error  is  defined  as  max  |c,|.  The  lz -error  is  the  appropriate 

measure  for  smooth  problems.  The  -error  is  more  appropriate  for  non-smooth  problems,  and 
should  be  the  primary  basis  for  comparison  in  these  tests.  Since  in  these  problems  the  largest  error 
occurs  at  discontinuities,  the  the  fj-error,  which  preferentially  measures  the  larger  errors,  will  be 
strongly  affected  by  very  small  changes  in  Cue  position  of  the  discontinuities.  The  l|-error,  on  the 
other  hand,  weights  all  errors  equally  and  is  a  better  measure  of  the  numerical  dispersion  away 
from  the  discontinuities.  The  ^-error  is  easily  contaminated  by  very  small  errors  in  positioning 
discontinuities  and  is  useful  only  for  continuous  solutions. 

As  a  first  sample  problem  we  take  a  variant  of  the  revolving  cone  problem  used  in  [6]  and  (2). 
The  problem  is: 

ui  -  yu,  +  xut  s=  0, 

with  the  initial  condition 


where 


«(*,  ¥,0) 


1  -  16r,  if  r  < 

0,  otherwise, 
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method 

ht 

*2 

error 

f-i 

error 

too 

error 

STEP 

Time  (seconds 
USINTG 

) 

BNDY 

coarse 

6.65  E-2 

2.19E-2 

5.98E-1 

3.73 

3.32 

0.33 

fine 

1/80 

1.38E-2 

3.74E-3 

1.22E-1 

209. 

204. 

4.0 

MR 

1.40E-2 

4.20E-3 

1.22E-1 

46.3 

35.0 

7.55 

Moving  MR 

1/80 

5.65E-3 

1.39E-3 

_ 

7.30E-2 

35.4 

25.8 

6.20 

TabU  1:  Results  for  first  test  problem.  Errors  are  computed 
on  the  coarse  grid.  Times  are  for  entire  program  (STEP), 
user  integrator  (USINTG),  and  boundary  handling  (BNDY). 


The  boundary  conditions  at  the  inflow  boundaries  are  zero.  The  solution  to  this  variable  coefficient 
problem  is  given  by  rigid  counter-clockwise  rotation  of  the  initial  data  about  the  origin  with  angular 
frequency  uj  =  1.  This  is  a  good  test  of  both  regular  and  moving  mesh  refinement  because  finite 
difference  approximations  to  this  problem  generate  large,  localized  errors  containing  high  frequency 
terms  (due  to  the  base  and  tip  of  the  cone).  For  moving  mesh  refinement,  there  is  the  additional 
feature  that  the  appropriate  velocity  varies  substantially  over  the  cone,  so  that  a  close  match  of 
the  velocity  of  the  solution  and  the  fine  grids  is  not  possible. 

In  our  test,  we  use  Lax-Wendroff  as  the  difference  approximation,  with  inflow  boundaries 
specified  as  u  =*  0  and  outflow  boundaries  specified  with  first  order  extrapolation.  A  substantial 
effort  was  made  to  identify  hotspots  in  the  code,  and  those  parts  of  the  cod<*  were  optimized.  This 
was  done  to  insure  that  these  tests  were  realistic  comparisons.  For  example,  the  execution  time 
of  the  Lax-Wendroff  scheme  was  improved  by  about  40%.  The  code  is  written  to  be  efficient  in 
all  modes  of  operation  (no  refinement,  unmoving  refinement,  or  moving  refinement),  and  is  fully 
instrumented. 

We  can  see  from  Table  1  that  the  MLUMR  solution  achieved  significantly  better  results  than 
LUMR  at  about  3/4  the  cost.  Both  LUMR  and  MLUMR  were  much  faster  than  a  uniform  grid. 
The  graphs  of  the  solutions  in  Figures  7-8  show  that  the  error  for  both  the  uniform  grid  and 
the  LUMR  calculation  is  due  primarily  to  numerical  dispersion  (the  high  frequency  waves  lag 
the  solution,  as  is  predicted  by  Figure  I),  while  the  error  in  the  MLUMR  solution  is  due  to  a 
combination  of  dispersive  error  (caused  by  the  variation  of  velocity  over  the  cone)  and  dissipation 
in  the  scheme.  We  can  see  from  Figure  2  that  part  of  the  efficiency  advantage  of  MLUMR  comes 
from  needing  less  refinement  (a  consequence  of  the  smaller  dispersive  error). 
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method 

hi 

h 

error 

*2 

error 

STEP 

Time  (seconds 
USINTG 

») 

BNDY 

coarse 

1/40 

1.84E-2 

4.78E-2 

.88 

.04 

fine 

1/80 

1.18E-2 

4.56E-2 

7.05 

.14 

MR 

1/80 

1.17E-2 

4.47E-2 

3.61 

.66 

Moving  MR 

1/40 

1.07E-2 

3.99E-2 

1.27 

.49 

Moving  MR 

1/80 

6.88E-3 

3.29E-2 

4.76 

2.97 

.97 

Table  3:  Results  for  second  test  problem.  Errors  are  com¬ 


puted  on  the  coarse  grid.  Times  are  for  entire  program 
(STEP),  user  integrator  (USINTG),  and  boundary  handling 
(BNDY). 

The  second  test  problem  is 


U|  +  uua  +  utiy  =  0 


on  the  unit  square,  with  0  <  t  <  0.6.  Note  that  under  the  rotation  x*  =  (x  +  y)/>/2,  y*  * 
(-z  +  y)/>/2,  the  equation  becomes  u(  +  y/2 uuy  =  0;  the  solutions  of  this  equation  are  well  known. 
Thu  is  a  very  difficult  problem  for  any  mesh  refinement  scheme  because  so  much  of  the  domain 
needs  to  be  refined. 

The  initial  data  for  our  second  test  problem  is 


u(x,y,t  =  0) 


if  0  <  z  <  j  and  0  <  y  <  |, 
-{,  if  |  <  z  <  1  and  5  <  y  <  1, 
7,  otherwise. 


The  solution  of  this  problem  consists  of  six  shocks,  two  moving  with  speed  3/8,  two  with  speed 
1/8,  and  two  stationary.  At  two  points  (the  problem  is  symmetric  about  z  ■  y),  three  shocks  with 
different  speeds  come  together.  Further,  since  this  is  a  nonlinear  problem,  even  a  perfect  choice  for 
the  speed  of  the  local  grids  will  not  completely  eliminate  the  error.  This  is  because  characteristics 
converge  on  the  shock,  and  thus  there  is  no  unique  local  velocity. 

For  this  problem,  the  grid  speeds  are  chosen  nearly  exactly;  the  speeds  were  set  by  the  location 
and  extent  of  the  grid.  If  the  grid  covered  several  shocks,  it  was  broken  up  into  several  overlapping 
pieces.  Because  the  grid  speeds  were  chosen  nearly  exactly,  the  errors  in  the  computation  come 
from  only  two  places:  the  errors  in  fitting  the  moving  grids  to  regions  where  multiple  shocks 
interact,  and  the  error  in  integrating  a  shock  even  in  a  co-moving  reference  frame.  Thus,  this 
example  is  complementary  to  the  first  example,  where  the  error  is  due  to  range  of  velocities  over 
the  refined  grid. 
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We  can  see  from  Table  2  that  even  without  refinement,  MLUMR  gives  a  much  better  solution 
at  a  fraction  of  the  cost  of  either  LUMR  or  uniform  grids.  In  this  example,  the  fi-norm  is  the 
most  representative,  as  an  error  of  one  mesh  point  at  a  she  can  make  the  4o*norm  0(1)  and 
have  significant  impact  on  the  fj-norm.  From  Figures  9-10  we  see  that  the  majority  of  the  error 
in  all  computations  is  at  the  shocks,  but  that  the  *ringing*  trailing  the  shocks  is  nearly  absent  in 
the  MLUMR  calculation.  The  highest  error  in  the  MLUMR  calculation  comes  at  the  intersections 
of  the  shocks,  where  the  MLUMR  method  can’t  take  advantage  of  any  local  frame  of  reference. 
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Coot's*  grid  solution  ot  T 


3.375 


Pin#  grid  solution  ot  T  s  3.375 


Flfort  7:  Solutions  for  revolving  cone  problem 


UJHft  solution  ot  T  :  0.600 


solution  ot  T  s  0.600 


Figure  9i  Solution  for  the  second  test  problem. 


Error  In  wrH  arid 


Error  In  ftno  grid 


Error  IU?  S. 6137-0011  Error  IUr  8.1310-0011 


Error  in  LUMH  oolution  Error  In  H.UHI  oolution 


Error  IUr  8.5516-0011  Error  (U=  7.6162-001 1 


Figure  IOi  Errors  in  the  solution*  in  Figure  9. 
16 


flfnra  11:  Grids  for  both  test  problems  at  the  final  time 
(7  »  3.375  for  problem  1,7  s  0.6  for  problem  2). 
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